/*******************************************************************************
 * Copyright (c) SCUT. 2022. All rights reserved.
 * Description: The realization of my_syr.
 * Author: CHEN Han
 * Create: 2022-06-26
 *******************************************************************************/

#include "my_syr.h"

void my_csyr(enum CBLAS_ORDER order,
             enum CBLAS_UPLO uplo,
             const int n,
             float complex Alpha,
             float complex* X,
             const int incx,
             float complex* A,
             const int lda)
             {
    if (order != CblasColMajor) {
        return;
    }

    if (n == 0 || Alpha == 0.0f) {
        return;
    }

    int i, info, ix, j, jx, kx;
    float complex temp;

    if (incx <= 0) {
        kx = 1 - (n - 1) * incx;
    } else if (incx != 1) {
        kx = 1;
    }

    if (uplo == CblasUpper) {
        if (incx == 1) {
            for (j = 1; j <= n; j++) {
                if (X[j - 1] != 0.0f) {
                    temp = Alpha * X[j - 1];
                    for (i = 1; i <= j; i++) {
                        A[i - 1 + (j - 1) * (lda)] += temp * X[i - 1];
                    }
                }
            }
        } else {
            jx = kx;
            for (j = 1; j <= n; j++) {
                if (X[jx - 1] != 0.0f) {
                    temp = Alpha * X[jx - 1];
                    ix = kx;
                    for (i = 1; i <= j; i++) {
                        A[i - 1 + (j - 1) * (lda)] += temp * X[ix - 1];
                        ix += incx;
                    }
                }
                jx += incx;
            }
        }
    } else {
        if (incx == 1) {
            for (j = 1; j <= n; j++) {
                if (X[j - 1] != 0.0) {
                    temp = Alpha * X[j - 1];
                    for (i = j; i <= n; i++) {
                        A[i - 1 + (j - 1) * (lda)] += temp * X[i - 1];
                    }
                }
            }
        } else {
            jx = kx;
            for (j = 1; j <= n; j++) {
                if (X[jx - 1] != 0.0f) {
                    temp = Alpha * X[jx - 1];
                    ix = kx;
                    for (i = j; i <= n; i++) {
                        A[i - 1 + (j - 1) * (lda)] += temp * X[ix - 1];
                        ix += incx;
                    }
                }
                jx += incx;
            }
        }
    }
}

void my_zsyr(enum CBLAS_ORDER order,
             enum CBLAS_UPLO uplo,
             const int n,
             double complex Alpha,
             double complex* X,
             const int incx,
             double complex* A,
             const int lda) {
    if (order != CblasColMajor) {
        exit(1);
    }

    if (n == 0 || Alpha == 0.0f) {
        return;
    }

    int i, info, ix, j, jx, kx;
    double complex temp;

    if (incx <= 0) {
        kx = 1 - (n - 1) * incx;
    } else if (incx != 1) {
        kx = 1;
    }

    if (uplo == CblasUpper) {
        if (incx == 1) {
            for (j = 1; j <= n; j++) {
                if (X[j - 1] != 0.0f) {
                    temp = Alpha * X[j - 1];
                    for (i = 1; i <= j; i++) {
                        A[i - 1 + (j - 1) * (lda)] += temp * X[i - 1];
                    }
                }
            }
        } else {
            jx = kx;
            for (j = 1; j <= n; j++) {
                if (X[jx - 1] != 0.0f) {
                    temp = Alpha * X[jx - 1];
                    ix = kx;
                    for (i = 1; i <= j; i++) {
                        A[i - 1 + (j - 1) * (lda)] += temp * X[ix - 1];
                        ix += incx;
                    }
                }
                jx += incx;
            }
        }
    } else {
        if (incx == 1) {
            for (j = 1; j <= n; j++) {
                if (X[j - 1] != 0.0f) {
                    temp = Alpha * X[j - 1];
                    for (i = j; i <= n; i++) {
                        A[i - 1 + (j - 1) * (lda)] += temp * X[i - 1];
                    }
                }
            }
        } else {
            jx = kx;
            for (j = 1; j <= n; j++) {
                if (X[jx - 1] != 0.0f) {
                    temp = Alpha * X[jx - 1];
                    ix = kx;
                    for (i = j; i <= n; i++) {
                        A[i - 1 + (j - 1) * (lda)] += temp * X[ix - 1];
                        ix += incx;
                    }
                }
                jx += incx;
            }
        }
    }
}